# 2016.11.08 重新整理cerna的结果，才发现有问题所以重新跑的结果
setwd("F:/cuffdiff")
# setwd("/home/zhang/BMK/test/cuffdiff")

library(openxlsx)

# 把结果读进来
cerna <- unique(read.table("cerna.txt", header = T))
cerna$q <- p.adjust(cerna$`pֵ值`,method = "BH" )

message <- read.table("~/index/name/allMessage")

# 整合rna1的信息
fir <- merge(cerna, message, by.x = "RNA_1", by.y = "V5")
fir <- unique(fir[,c(1,9,10,11,12,14,15, 2:8)])
colnames(fir) <- c("RNA_1", "chr_1", "start_1", "end_1", "strand_1", "type_1", "gene_name_1", "RNA_2", "N", "K", "n", "c", "p.value", "q.value")
# 整合RNA2的信息
sec <- merge(fir, message, by.x = "RNA_2", by.y = "V5")
sec <- unique(sec[,c(2:8, 1, 15:18, 20,21, 9:14)])
colnames(sec)[9:14] <- c("chr_2", "start_2", "end_2", "strand_2", "type_2", "gene_name_2")
# 抽取有差异的
diff <- unique(sec[sec$q.value < 0.05,])
# 写出来
write.table(sec, "cerna_allmessage.txt", row.names = F, col.names = T, quote = F,sep = "\t")
write.csv(sec, "cerna.csv", row.names = F)
write.csv(diff, "cerna_diff.csv", row.names = F)
write.xlsx(diff[diff$q.value < 0.01,], "cerna_superdiff.xlsx", rowNames = F)


diff <- read.csv("cerna结果/cerna_diff.csv")
# 提取所有后边要用到的miranda
miranda <- read.table("all.miranda")

change.type <- function(input){
  input <- input[3]
  if(input == "protein_coding"){
    input = "mrna"
  }else{
    input = "lnc"
  }
  return(input)
}


# 提取rna的数据
rna <- diff[,c(2:4,7,6,5,1)]
colnames(rna) <- colnames(diff)[c(9:11,14,13,12,8)]
rna <- rbind(rna, diff[,c(9:11,14,13,12,8)])
rna <- unique(rna)

test <- merge(miranda, rna, by.x = "V2", by.y = "RNA_2")
test <- unique(test[,c(2,1,7)])
test[,3] <- apply(test, 1, change.type)
# 这是q < 0.05的
# 生成mrna列表和lnc的bed格式，用于后来的lnc功能注释
mrna <- unique(subset(rna, type_2 == "protein_coding"))
lnc <- unique(subset(rna, type_2 != "protein_coding"))


write.table(mrna$gene_name_2, "mrna.cerna", quote = F, row.names = F, col.names = F)
write.table(lnc[,c(1:3,7,2,6)], "bed文件/lnc_cerna.bed", row.names = F, col.names = F, quote = F, sep = "\t")
write.table(test, "diff.miranda", row.names = F, col.names = F, quote = F, sep = "\t")
# 生成superdiff的mrna的列表和lnc的bed格式
diff <- diff[diff$q.value < 0.01,]
rna <- diff[,c(2:4,7,6,5,1)]
colnames(rna) <- colnames(diff)[c(9:11,14,13,12,8)]
rna <- rbind(rna, diff[,c(9:11,14,13,12,8)])
rna <- unique(rna)

test <- merge(miranda, rna, by.x = "V2", by.y = "RNA_2")
test <- unique(test[,c(2,1,7)])
test[,3] <- apply(test, 1, change.type)

mrna <- unique(subset(rna, type_2 == "protein_coding"))
lnc <- unique(subset(rna, type_2 != "protein_coding"))

write.table(mrna$gene_name_2, "mrna_superdiff.cerna", quote = F, row.names = F, col.names = F)
write.table(lnc[,c(1:3,7,2,6)], "bed文件/lnc_superdiff_cerna.bed", row.names = F, col.names = F, quote = F, sep = "\t")
write.table(test, "superdiff.miranda", row.names = F, col.names = F, quote = F, sep = "\t")

